*************************************************************
*** Monte Carlo experiments (SLX DGP) for the WWW JoP Project
***
*** Created: 9-19-14
*** Modified: 9-5-19
***
*************************************************************

clear 
clear matrix
clear mata
set matsize 1000
set maxvar 32000

version 14

do "Programs\SAR Programs.do"

*************************************************************
*** Spatial heterogeneity (varying coefficients for thetas)
*************************************************************

*************************************************************
*** Set up the parameters of the MC experiments
*************************************************************
matrix beta = (1 \ 1)		/* Parameters for beta_x1, and beta_x2 */
matrix I = I(200)			/* Identity matrix with N observations */

local s = 2					/* Standard deviation of the errors */
local MCsims = 1000			/* Number of simulations */
local QIsims = 1000
local obs = colsof(I)		/* Number of observations */
local seed = 648			/* Set the seed */

set seed `seed'
set obs `obs'
if `obs' > `QIsims' {
	set matsize `obs'
}
else {
	set matsize `QIsims'
}
*************************************************************

spatwmat using "Data\W\W_`obs'rs.dta", name(W)

*** Coefficients for theta: WX_1 (row 1) and WX_2 (row 2)
mat Theta = (-0.8, -0.8, -0.8, 0.4,  0.4,  0.4, -0.2, -0.2, -0.2 \ /*
*/			 -0.8, 0.01,  0.4, 0.4, 0.01, -0.8, -0.2, 0.01,  0.4)											
local end = colsof(Theta)

mat W2 = W*W

mata: I = st_matrix("I")
mata: W = st_matrix("W")
mata: W2 = st_matrix("W2")

tempvar id
gen `id' = _n

preserve
	use "Data\W\W_`obs'rs.dta", clear
	gen id = _n
	spmat dta Wspat V*, id(id) replace
restore

*************************************************************
*** Calculate the True Effects
*************************************************************

tempname slxtrue
postfile `slxtrue' b_x1 b_x2 theta_1 theta_2  /* 
*/	true_teffect_x1 true_deffect_x1 true_ieffect_x1 true_feffect_x1 true_heffect_x1 /*
*/	true_teffect_x2 true_deffect_x2 true_ieffect_x2 true_feffect_x2 true_heffect_x2 /*
*/	true_deffect_o0_x1 /*
*/	true_deffect_o0_x2 /*
*/	true_teffect_o1_x1 true_deffect_o1_x1 true_ieffect_o1_x1 /*
*/	true_teffect_o1_x2 true_deffect_o1_x2 true_ieffect_o1_x2 /*
*/	true_teffect_o2_x1 true_deffect_o2_x1 true_ieffect_o2_x1 /*
*/	true_teffect_o2_x2 true_deffect_o2_x2 true_ieffect_o2_x2 /*
*/	true_teffect_o3_x1 true_deffect_o3_x1 true_ieffect_o3_x1 /*
*/	true_teffect_o3_x2 true_deffect_o3_x2 true_ieffect_o3_x2 /*
*/	using "Generate\SLX\SLX--True Effects.dta", replace every(20)

qui foreach c of numlist 1(1)`end' {
	foreach z of numlist 1 2 {
		scalar beta_x`z' = beta[`z',1]
		scalar theta_x`z' = Theta[`z', `c']
		
		matrix true = [beta[1,1], Theta[`z', `c']]
		mata: true = st_matrix("true")
		mata: calceffslx1(true, W, 1)
	
		foreach e in t d i {	
			scalar true_`e'effect_x`z' = r(`e'effect)
			scalar true_`e'effect_o1_x`z' = r(`e'effect_o1)
			scalar true_`e'effect_o2_x`z' = 0
			scalar true_`e'effect_o3_x`z' = 0			
		}

		foreach e in feffect heffect {
			scalar true_`e'_x`z' = 0
		}

		scalar true_deffect_o0_x`z' = r(deffect_o0)
	}
	
	post `slxtrue' (beta_x1) (beta_x2) (theta_x1) (theta_x2) /* 
	*/	(true_teffect_x1) (true_deffect_x1) (true_ieffect_x1) (true_feffect_x1) (true_heffect_x1) /*  
	*/	(true_teffect_x2) (true_deffect_x2) (true_ieffect_x2) (true_feffect_x2) (true_heffect_x2) /* 
	*/	(true_deffect_o0_x1) /*
	*/	(true_deffect_o0_x2) /*
	*/	(true_teffect_o1_x1) (true_deffect_o1_x1) (true_ieffect_o1_x1) /*
	*/	(true_teffect_o1_x2) (true_deffect_o1_x2) (true_ieffect_o1_x2) /*
	*/	(true_teffect_o2_x1) (true_deffect_o2_x1) (true_ieffect_o2_x1) /*
	*/	(true_teffect_o2_x2) (true_deffect_o2_x2) (true_ieffect_o2_x2) /*
	*/	(true_teffect_o3_x1) (true_deffect_o3_x1) (true_ieffect_o3_x1) /*
	*/	(true_teffect_o3_x2) (true_deffect_o3_x2) (true_ieffect_o3_x2) 	
}
	
postclose `slxtrue'	

preserve
	use "Generate\SLX\SLX--True Effects.dta", clear
	sort theta_1 theta_2
	save "Generate\SLX\SLX--True Effects.dta", replace
restore


*************************************************************
*** Monte Carlo Experiments
*************************************************************
tempname slx
postfile `slx' sim str30 time theta_1 theta_2 /*
*/	b_x1_m1 se_x1_m1 b_x2_m1 se_x2_m1 teffect_x1_m1 teffect_x1_se_m1 deffect_x1_m1 deffect_x1_se_m1 deffect_o0_x1_m1 deffect_o0_x1_se_m1 teffect_x2_m1 teffect_x2_se_m1 deffect_x2_m1 deffect_x2_se_m1 deffect_o0_x2_m1 deffect_o0_x2_se_m1 ar2_m1 vif_m1 /*
*/	b_x1_m2 se_x1_m2 b_x2_m2 se_x2_m2 theta_1_m2 se_theta_1_m2 theta_2_m2 se_theta_2_m2 p_theta_1_m2 p_theta_2_m2 p_theta_diff_m2 p_theta_joint_m2 ar2_m2 vif_m2 /*
*/	teffect_x1_m2 teffect_x1_se_m2 deffect_x1_m2 deffect_x1_se_m2 ieffect_x1_m2 ieffect_x1_se_m2 deffect_o0_x1_m2 deffect_o0_x1_se_m2 teffect_o1_x1_m2 teffect_o1_x1_se_m2 deffect_o1_x1_m2 deffect_o1_x1_se_m2 ieffect_o1_x1_m2 ieffect_o1_x1_se_m2 /*
*/	teffect_x2_m2 teffect_x2_se_m2 deffect_x2_m2 deffect_x2_se_m2 ieffect_x2_m2 ieffect_x2_se_m2 deffect_o0_x2_m2 deffect_o0_x2_se_m2 teffect_o1_x2_m2 teffect_o1_x2_se_m2 deffect_o1_x2_m2 deffect_o1_x2_se_m2 ieffect_o1_x2_m2 ieffect_o1_x2_se_m2 /*
*/	b_x1_m3 se_x1_m3 b_x2_m3 se_x2_m3 theta_1_m3 se_theta_1_m3 theta_2_m3 se_theta_2_m3 p_theta_1_m3 p_theta_2_m3 p_theta_diff_m3 ar2_m3 vif_m3 /*
*/	teffect_x1_m3 teffect_x1_se_m3 deffect_x1_m3 deffect_x1_se_m3 ieffect_x1_m3 ieffect_x1_se_m3 feffect_x1_m3 feffect_x1_se_m3 heffect_x1_m3 heffect_x1_se_m3 deffect_o0_x1_m3 deffect_o0_x1_se_m3 teffect_o1_x1_m3 teffect_o1_x1_se_m3 deffect_o1_x1_m3 deffect_o1_x1_se_m3 ieffect_o1_x1_m3 ieffect_o1_x1_se_m3 teffect_o2_x1_m3 teffect_o2_x1_se_m3 deffect_o2_x1_m3 deffect_o2_x1_se_m3 ieffect_o2_x1_m3 ieffect_o2_x1_se_m3 /*
*/	teffect_x2_m3 teffect_x2_se_m3 deffect_x2_m3 deffect_x2_se_m3 ieffect_x2_m3 ieffect_x2_se_m3 feffect_x2_m3 feffect_x2_se_m3 heffect_x2_m3 heffect_x2_se_m3 deffect_o0_x2_m3 deffect_o0_x2_se_m3 teffect_o1_x2_m3 teffect_o1_x2_se_m3 deffect_o1_x2_m3 deffect_o1_x2_se_m3 ieffect_o1_x2_m3 ieffect_o1_x2_se_m3 teffect_o2_x2_m3 teffect_o2_x2_se_m3 deffect_o2_x2_m3 deffect_o2_x2_se_m3 ieffect_o2_x2_m3 ieffect_o2_x2_se_m3 /*
*/	b_x1_m4 se_x1_m4 b_x2_m4 se_x2_m4 rho_m4 rho_se_m4 rmse_m4 /*
*/	teffect_x1_m4 teffect_x1_se_m4 deffect_x1_m4 deffect_x1_se_m4 ieffect_x1_m4 ieffect_x1_se_m4 feffect_x1_m4 feffect_x1_se_m4 heffect_x1_m4 heffect_x1_se_m4 deffect_o0_x1_m4 deffect_o0_x1_se_m4 /*
*/	teffect_o1_x1_m4 teffect_o1_x1_se_m4 deffect_o1_x1_m4 deffect_o1_x1_se_m4 ieffect_o1_x1_m4 ieffect_o1_x1_se_m4 teffect_o2_x1_m4 teffect_o2_x1_se_m4 deffect_o2_x1_m4 deffect_o2_x1_se_m4 ieffect_o2_x1_m4 ieffect_o2_x1_se_m4 teffect_o3_x1_m4 teffect_o3_x1_se_m4 deffect_o3_x1_m4 deffect_o3_x1_se_m4 ieffect_o3_x1_m4 ieffect_o3_x1_se_m4 /*
*/	teffect_x2_m4 teffect_x2_se_m4 deffect_x2_m4 deffect_x2_se_m4 ieffect_x2_m4 ieffect_x2_se_m4 feffect_x2_m4 feffect_x2_se_m4 heffect_x2_m4 heffect_x2_se_m4 deffect_o0_x2_m4 deffect_o0_x2_se_m4 /* 
*/	teffect_o1_x2_m4 teffect_o1_x2_se_m4 deffect_o1_x2_m4 deffect_o1_x2_se_m4 ieffect_o1_x2_m4 ieffect_o1_x2_se_m4 teffect_o2_x2_m4 teffect_o2_x2_se_m4 deffect_o2_x2_m4 deffect_o2_x2_se_m4 ieffect_o2_x2_m4 ieffect_o2_x2_se_m4 teffect_o3_x2_m4 teffect_o3_x2_se_m4 deffect_o3_x2_m4 deffect_o3_x2_se_m4 ieffect_o3_x2_m4 ieffect_o3_x2_se_m4 /*
*/	using "Generate\SLX\SLX.dta", replace every(10)

local stfull = "$S_TIME"
local operfull = "SLX DGP: `MCsims' simulations of `end' scenarios"

local b_1 = beta[1,1]
local b_2 = beta[2,1]
		
foreach s of numlist 1(1)`end' {
	local theta_1 = Theta[1,`s']
	local theta_2 = Theta[2,`s']

	qui foreach i of numlist 1(1)`MCsims' {
		local st1 = "$S_TIME"
		local oper1 = "Generate the variables and estimate SLX models"
		
		drawnorm r, n(`obs') sd(`s')
		generate x1 = -10+(10 - -10)*uniform()
		generate x2 = -5+(5 - -5)*uniform()

		*** Generate the two spatial-lags
		foreach m of numlist 1 2 {
			mkmat x`m'
			
			mat slag`m' = W * x`m'
			svmat slag`m'
			rename slag`m'1 slag`m'
			
			mat slag`m'_2 = W2 * x`m'
			svmat slag`m'_2
			rename slag`m'_21 slag`m'_2			
		}

		*=========================================================
		*=========================================================
		*** Generate the SLX data-generating process
		gen y = (`b_1' * x1) + (`b_2' * x2) + (`theta_1' * slag1) + (`theta_2' * slag2) + r

		*=========================================================
		*=========================================================
		
		if `obs' > `QIsims' {
			set obs `obs'
}
		else {
			set obs `QIsims'
		}	
		
		* Model 1: under-specified SLX
		local m = 1
		
		qui reg y x1 x2 
		foreach x of numlist 1 2{
			scalar b_x`x'_m`m' = _b[x`x']
			scalar se_x`x'_m`m' = _se[x`x']
			
			foreach v in teffect deffect deffect_o0 {
				scalar `v'_x`x'_m`m' = b_x`x'_m`m'
				scalar `v'_x`x'_se_m`m' = se_x`x'_m`m'
			}
		}

		scalar ar2_m`m' = e(r2_a)
		vif 
		scalar vif_m`m' = r(vif_1)

		* Model 2: standard SLX (correct model)
		local m = 2
		
		qui reg y x1 x2 slag1 slag2
		mat V = e(V)

		foreach k of numlist 1 2{
			scalar b_x`k'_m`m' = _b[x`k']
			scalar se_x`k'_m`m' = _se[x`k']
		
			scalar theta_x`k'_m`m' = _b[slag`k']
			scalar se_theta_x`k'_m`m' = _se[slag`k']
		
			test slag`k' = 0
			scalar p_theta_`k'_m`m' = r(p)
		}

		test slag1 = slag2 
		scalar p_theta_diff_m`m' = r(p)

		test slag1 = slag2 = 0
		scalar p_theta_joint_m`m' = r(p)
	
		scalar ar2_m`m' = e(r2_a)
		vif
		scalar vif_m`m' = r(vif_1)		
		
		matrix coeff = [_b[x1], _b[x2], _b[slag1], _b[slag2]]
		matrix VCVM = [V[1..4,1..4]]
		
		preserve	
			clear
			tempvar beta1_draw beta2_draw theta1_draw theta2_draw 
			corr2data `beta1_draw' `beta2_draw' `theta1_draw' `theta2_draw', n(`QIsims') means(coeff) cov(VCVM)
			mkmat `beta1_draw' `beta2_draw' `theta1_draw' `theta2_draw', matrix(draws_m`m')	
		restore		
		
		mat draws1_m`m' = [draws_m`m'[1...,1], draws_m`m'[1...,3]]
		mata: draws1_m`m' = st_matrix("draws1_m`m'")
		mata: calceffslx1(draws1_m`m', W, `QIsims')

		local x = 1
		scalar deffect_o0_x`x'_m`m' = r(deffect_o0)
		scalar deffect_o0_x`x'_se_m`m' = r(deffect_se_o0)
		
		foreach e in t d i {	
			scalar `e'effect_x`x'_m`m' = r(`e'effect)
			scalar `e'effect_x`x'_se_m`m' = r(`e'effect_se)
			
			scalar `e'effect_o1_x`x'_m`m' = r(`e'effect_o1)
			scalar `e'effect_o1_x`x'_se_m`m' = r(`e'effect_se_o1)
		}		
	
		mat draws2_m`m' = [draws_m`m'[1...,2], draws_m`m'[1...,4]]
		mata: draws2_m`m' = st_matrix("draws2_m`m'")
		mata: calceffslx1(draws2_m`m', W, `QIsims')

		local x = 2
		scalar deffect_o0_x`x'_m`m' = r(deffect_o0)
		scalar deffect_o0_x`x'_se_m`m' = r(deffect_se_o0)
		
		foreach e in t d i {	
			scalar `e'effect_x`x'_m`m' = r(`e'effect)
			scalar `e'effect_x`x'_se_m`m' = r(`e'effect_se)
			
			scalar `e'effect_o1_x`x'_m`m' = r(`e'effect_o1)
			scalar `e'effect_o1_x`x'_se_m`m' = r(`e'effect_se_o1)
		}		
		
		* Model 3: overspecified complete model
		local m = 3
		
		qui reg y x1 x2 slag1 slag2 slag1_2 slag2_2
		mat b = e(b)
		mat V = e(V)
		
		matrix coeff = b[1...,1..6]
		matrix VCVM = [V[1..6,1..6]]
		
		preserve	
			clear
			tempvar beta1_draw beta2_draw theta1_1_draw theta2_1_draw theta1_2_draw theta2_2_draw
			corr2data `beta1_draw' `beta2_draw' `theta1_1_draw' `theta2_1_draw' `theta1_2_draw' `theta2_2_draw', n(`QIsims') means(coeff) cov(VCVM)
			mkmat `beta1_draw' `beta2_draw' `theta1_1_draw' `theta2_1_draw' `theta1_2_draw' `theta2_2_draw', matrix(draws_m`m')	
		restore		
		
		mat draws1_m`m' = [draws_m`m'[1...,1], draws_m`m'[1...,3], draws_m`m'[1...,5]]
		mata: draws1_m`m' = st_matrix("draws1_m`m'")
		mata: calceffslx2(draws1_m`m', W, W2, `QIsims')

		local x = 1
		foreach e in t d i {	
			foreach o in 1 2 {
				scalar `e'effect_o`o'_x`x'_m`m' = r(`e'effect_o`o')
				scalar `e'effect_o`o'_x`x'_se_m`m' = r(`e'effect_se_o`o')
			}
		}

		foreach e in t d i f h {
			scalar `e'effect_x`x'_m`m' = r(`e'effect)
			scalar `e'effect_x`x'_se_m`m' = r(`e'effect_se)		
		}

		scalar deffect_o0_x`x'_m`m' = r(deffect_o0)
		scalar deffect_o0_x`x'_se_m`m' = r(deffect_se_o0)
		
		mat draws2_m`m' = [draws_m`m'[1...,2], draws_m`m'[1...,4], draws_m`m'[1...,6]]
		mata: draws2_m`m' = st_matrix("draws2_m`m'")
		mata: calceffslx2(draws2_m`m', W, W2, `QIsims')

		local x = 2
		foreach e in t d i {	
			foreach o in 1 2 {
				scalar `e'effect_o`o'_x`x'_m`m' = r(`e'effect_o`o')
				scalar `e'effect_o`o'_x`x'_se_m`m' = r(`e'effect_se_o`o')
			}
		}

		foreach e in t d i f h {
			scalar `e'effect_x`x'_m`m' = r(`e'effect)
			scalar `e'effect_x`x'_se_m`m' = r(`e'effect_se)		
		}

		scalar deffect_o0_x`x'_m`m' = r(deffect_o0)
		scalar deffect_o0_x`x'_se_m`m' = r(deffect_se_o0)
		
		foreach x of numlist 1 2{
			scalar b_x`x'_m`m' = _b[x`x']
			scalar se_x`x'_m`m' = _se[x`x']
		
			scalar theta_x`x'_m`m' = _b[slag`x']
			scalar se_theta_x`x'_m`m' = _se[slag`x']
		
			test slag`x' = 0
			scalar p_theta_`x'_m`m' = r(p)
		}

		test slag1 = slag2 
		scalar p_theta_diff_m`m' = r(p)
	
		scalar ar2_m`m' = e(r2_a)
		vif
		scalar vif_m`m' = r(vif_1)
	
		* Model 4: SAR Model
		local m = 4
		
		spreg ml y x1 x2 in 1/`obs', id(`id') dlmat(Wspat, eig)
		mat b = e(b)
		mat V = e(V)
				
		matrix coeff = [b[1,4], b[1,1]]
		matrix VCVM = [V[4,4], V[4,1] \ V[1,4], V[1,1]]
		
		preserve
			clear
			tempvar rho_draw beta_draw
			corr2data `rho_draw' `beta_draw', n(`QIsims') means(coeff) cov(VCVM)
			mkmat `rho_draw' `beta_draw', matrix(draws1_sar)
		restore
		
		mata: draws1_sar = st_matrix("draws1_sar")
		mata: calceffsar(draws1_sar, I, W, `QIsims')

		local x = 1
		foreach e in t d i {
			scalar `e'effect_x`x'_m4 = r(`e'effect)
			scalar `e'effect_x`x'_se_m4 = r(`e'effect_se)		
		
			foreach o of numlist 1(1)3 {
				scalar `e'effect_o`o'_x`x'_m4 = r(`e'effect`o')
				scalar `e'effect_o`o'_x`x'_se_m4 = r(`e'effect`o'_se)
			}
		}
		
		foreach e in f h {
			scalar `e'effect_x`x'_m4 = r(`e'effect)
			scalar `e'effect_x`x'_se_m4 = r(`e'effect_se)
		}
		
		scalar deffect_o0_x`x'_m4 = r(deffect0)
		scalar deffect_o0_x`x'_se_m4 = r(deffect0_se)
		
		matrix coeff = [b[1,4], b[1,2]]
		matrix VCVM = [V[4,4], V[4,2] \ V[2,4], V[2,2]]
		
		preserve
			clear
			tempvar rho_draw beta_draw
			corr2data `rho_draw' `beta_draw', n(`QIsims') means(coeff) cov(VCVM)
			mkmat `rho_draw' `beta_draw', matrix(draws2_sar)
		restore
		
		mata: draws2_sar = st_matrix("draws2_sar")
		mata: calceffsar(draws2_sar, I, W, `QIsims')
		
		local x = 2
		foreach e in t d i {
			scalar `e'effect_x`x'_m4 = r(`e'effect)
			scalar `e'effect_x`x'_se_m4 = r(`e'effect_se)		
		
			foreach o of numlist 1(1)3 {
				scalar `e'effect_o`o'_x`x'_m4 = r(`e'effect`o')
				scalar `e'effect_o`o'_x`x'_se_m4 = r(`e'effect`o'_se)
			}
		}
		
		foreach e in f h {
			scalar `e'effect_x`x'_m4 = r(`e'effect)
			scalar `e'effect_x`x'_se_m4 = r(`e'effect_se)
		}
		
		scalar deffect_o0_x`x'_m4 = r(deffect0)
		scalar deffect_o0_x`x'_se_m4 = r(deffect0_se)
		
		foreach x of numlist 1 2 {
			scalar b_x`x'_m`m' = _b[y:x`x']
			scalar se_x`x'_m`m' = _se[y:x`x']		
		
		}
		
		scalar rho_m`m' = _b[lambda:_cons]
		scalar rho_se_m`m' = _se[lambda:_cons]

		spreg_fit
		scalar rmse_m`m' = round(r(rmse), .001)		

		if mod(`i',5) == 0 {
			nois display "." _c
			if mod(`i',100) == 0 {
				nois display "$S_DATE : `c(current_time)': Simulation #`i'"
			}
		}		

		elapse "`st1'" "`oper1'"
		
		post `slx' (`i') ("$S_elap") (`theta_1') (`theta_2') /*
		*/	(b_x1_m1) (se_x1_m1) (b_x2_m1) (se_x2_m1) (teffect_x1_m1) (teffect_x1_se_m1) (deffect_x1_m1) (deffect_x1_se_m1) (deffect_o0_x1_m1) (deffect_o0_x1_se_m1) (teffect_x2_m1) (teffect_x2_se_m1) (deffect_x2_m1) (deffect_x2_se_m1) (deffect_o0_x2_m1) (deffect_o0_x2_se_m1) (ar2_m1) (vif_m1) /*
		*/	(b_x1_m2) (se_x1_m2) (b_x2_m2) (se_x2_m2) (theta_x1_m2) (se_theta_x1_m2) (theta_x2_m2) (se_theta_x2_m2) (p_theta_1_m2) (p_theta_2_m2) (p_theta_diff_m2) (p_theta_joint_m2) (ar2_m2) (vif_m2) /*
		*/	(teffect_x1_m2) (teffect_x1_se_m2) (deffect_x1_m2) (deffect_x1_se_m2) (ieffect_x1_m2) (ieffect_x1_se_m2) (deffect_o0_x1_m2) (deffect_o0_x1_se_m2) (teffect_o1_x1_m2) (teffect_o1_x1_se_m2) (deffect_o1_x1_m2) (deffect_o1_x1_se_m2) (ieffect_o1_x1_m2) (ieffect_o1_x1_se_m2) /*
		*/	(teffect_x2_m2) (teffect_x2_se_m2) (deffect_x2_m2) (deffect_x2_se_m2) (ieffect_x2_m2) (ieffect_x2_se_m2) (deffect_o0_x2_m2) (deffect_o0_x2_se_m2) (teffect_o1_x2_m2) (teffect_o1_x2_se_m2) (deffect_o1_x2_m2) (deffect_o1_x2_se_m2) (ieffect_o1_x2_m2) (ieffect_o1_x2_se_m2) /*		
		*/	(b_x1_m3) (se_x1_m3) (b_x2_m3) (se_x2_m3) (theta_x1_m3) (se_theta_x1_m3) (theta_x2_m3) (se_theta_x2_m3) (p_theta_1_m3) (p_theta_2_m3) (p_theta_diff_m3) (ar2_m3) (vif_m3) /*
		*/	(teffect_x1_m3) (teffect_x1_se_m3) (deffect_x1_m3) (deffect_x1_se_m3) (ieffect_x1_m3) (ieffect_x1_se_m3) (feffect_x1_m3) (feffect_x1_se_m3) (heffect_x1_m3) (heffect_x1_se_m3) (deffect_o0_x1_m3) (deffect_o0_x1_se_m3) (teffect_o1_x1_m3) (teffect_o1_x1_se_m3) (deffect_o1_x1_m3) (deffect_o1_x1_se_m3) (ieffect_o1_x1_m3) (ieffect_o1_x1_se_m3) (teffect_o2_x1_m3) (teffect_o2_x1_se_m3) (deffect_o2_x1_m3) (deffect_o2_x1_se_m3) (ieffect_o2_x1_m3) (ieffect_o2_x1_se_m3) /*
		*/	(teffect_x2_m3) (teffect_x2_se_m3) (deffect_x2_m3) (deffect_x2_se_m3) (ieffect_x2_m3) (ieffect_x2_se_m3) (feffect_x2_m3) (feffect_x2_se_m3) (heffect_x2_m3) (heffect_x2_se_m3) (deffect_o0_x2_m3) (deffect_o0_x2_se_m3) (teffect_o1_x2_m3) (teffect_o1_x2_se_m3) (deffect_o1_x2_m3) (deffect_o1_x2_se_m3) (ieffect_o1_x2_m3) (ieffect_o1_x2_se_m3) (teffect_o2_x2_m3) (teffect_o2_x2_se_m3) (deffect_o2_x2_m3) (deffect_o2_x2_se_m3) (ieffect_o2_x2_m3) (ieffect_o2_x2_se_m3) /*
		*/	(b_x1_m4) (se_x1_m4) (b_x2_m4) (se_x2_m4) (rho_m4) (rho_se_m4) (rmse_m4) /*
		*/	(teffect_x1_m4) (teffect_x1_se_m4) (deffect_x1_m4) (deffect_x1_se_m4) (ieffect_x1_m4) (ieffect_x1_se_m4) (feffect_x1_m4) (feffect_x1_se_m4) (heffect_x1_m4) (heffect_x1_se_m4) (deffect_o0_x1_m4) (deffect_o0_x1_se_m4) /*
		*/	(teffect_o1_x1_m4) (teffect_o1_x1_se_m4) (deffect_o1_x1_m4) (deffect_o1_x1_se_m4) (ieffect_o1_x1_m4) (ieffect_o1_x1_se_m4) (teffect_o2_x1_m4) (teffect_o2_x1_se_m4) (deffect_o2_x1_m4) (deffect_o2_x1_se_m4) (ieffect_o2_x1_m4) (ieffect_o2_x1_se_m4) (teffect_o3_x1_m4) (teffect_o3_x1_se_m4) (deffect_o3_x1_m4) (deffect_o3_x1_se_m4) (ieffect_o3_x1_m4) (ieffect_o3_x1_se_m4) /*
		*/	(teffect_x2_m4) (teffect_x2_se_m4) (deffect_x2_m4) (deffect_x2_se_m4) (ieffect_x2_m4) (ieffect_x2_se_m4) (feffect_x2_m4) (feffect_x2_se_m4) (heffect_x2_m4) (heffect_x2_se_m4) (deffect_o0_x2_m4) (deffect_o0_x2_se_m4) /*
		*/	(teffect_o1_x2_m4) (teffect_o1_x2_se_m4) (deffect_o1_x2_m4) (deffect_o1_x2_se_m4) (ieffect_o1_x2_m4) (ieffect_o1_x2_se_m4) (teffect_o2_x2_m4) (teffect_o2_x2_se_m4) (deffect_o2_x2_m4) (deffect_o2_x2_se_m4) (ieffect_o2_x2_m4) (ieffect_o2_x2_se_m4) (teffect_o3_x2_m4) (teffect_o3_x2_se_m4) (deffect_o3_x2_m4) (deffect_o3_x2_se_m4) (ieffect_o3_x2_m4) (ieffect_o3_x2_se_m4)
		keep in 1/`obs'	
		drop r x1 x2 y slag* 
		mata: mata drop draws1_m2 draws2_m2 draws1_m3 draws2_m3 draws1_sar draws2_sar
	}
}

postclose `slx'
elapse "`stfull'" "`operfull'"
